#在服务器上运行代码
setwd("/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/20210324")
library(openxlsx)
filea=read.csv("all.FDR.sig.at.least.one.add.direction.same.diff.csv",head=T)
filea$id=paste(filea$Chr,filea$Start,sep = ":")
filea1=filea[filea$FDR.sig>1,]
file=read.table("53K.add.GWAS.eQTL.DEG.motif.for.analysis.txt",header=T,sep="\t")
file=file[!duplicated(file$id),]	#对TF去重
file1=file[file$pattern.not.rm.dupl.num.DC>1,]
file2=file1[file1$BF_in_DC>1,]
file3=file1[file1$BF_in_DC>10,]

asm2=read.table("220520ASMs_anno.hg19_multianno.csv",head=T,sep=",")	#220K
asm2$unitID=paste(asm2$Chr,asm2$Start,asm2$Ref,asm2$Alt,sep=":")
asm=read.csv("869727.all.snp.vaf.up.down.20210321.csv",head=T)
asm=asm[asm$unitID %in% asm2$unitID,]
asm=asm[as.character(asm$unitID) %in% as.character(filea$unitID),]
#asmdata=filea1[filea$unitID %in% asm$unitID,]
asm1=tidyr::separate(asm,unitID,into=c("chr","pos","ref","alt"),sep=":")
asm1$id=paste(asm1$chr,asm1$pos,sep=":")

###control
con=read.table("all.heter.snp",head=T,sep="\t")
con=con[!con$id %in% as.character(filea$id),]		###remove the ASH
conid=as.character(con$id)
conid.random=sample(conid,150000,replace = FALSE)

#GTEx eQTL

eqtl=fread("/mnt/md1200/5/zhaocunyou/wangzhongju/eQTL_analysis_by_wangzhongju/eQTL_data/GTEx_v6p_eQTL_for_lqy.txt",header=T,sep="\t")
eqtl1=tidyr::separate(eqtl,SNP,into=c("Chr","Pos","Alt","Ref","Other"),sep="_")
eqtl1$Chr=paste0("chr",eqtl1$Chr)
eqtl1$id=paste(eqtl1$Chr,eqtl1$Pos,sep=":")
GTEx =data.frame(id=eqtl1$id,GTx.p.value=eqtl1$P,GTx.eQTL.tissue=eqtl1$tissue)

###							分不同组织进行富集分析，与LIBD.eQTL数据形成对比，相互佐证。
str1 =as.character(unique(GTEx$GTx.eQTL.tissue))

for(randomid in 1:3){
conid.random=sample(conid,150000,replace = FALSE)
con1=con[con$id %in% conid.random,]
con1=merge(con1,anno,by="id",all.x=T)

for(tis in str1){
gtx=GTEx[GTEx$GTx.eQTL.tissue==tis,]

result=data.frame(matrix(,nrow = 1,ncol =4))
names(result)=c("Term","overlap.num","OR","P.value")
result =result[-1,]


c=length(intersect(as.character(con1$id),as.character(gtx$id)))
d=150000-c
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
rt[1,]=c("control",c,NA,NA)
result=rbind(result,rt)

a=length(intersect(as.character(filea$id),as.character(gtx$id)))
b=dim(filea)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("117012",a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(filea1$id),as.character(gtx$id)))
b=dim(filea1)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("61725",a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(file$id),as.character(gtx$id)))
b=dim(file)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("53425",a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(file1$id),as.character(gtx$id)))
b=dim(file1)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("8544",a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(file2$id),as.character(gtx$id)))
b=dim(file2)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("807",a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(file3$id),as.character(gtx$id)))
b=dim(file3)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("200",a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(asm1$id),as.character(gtx$id)))
b=dim(asm1)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("13649.117KASH.22OKASM",a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

rtn =paste("./GTEx.result/GTEx.eQTL.",gsub("_Analysis.v6p","",tis),".enrich.random",randomid,".csv",sep="")
write.csv(result,rtn,quote=F,row.names=F)
}

}

###对总的GTEx eQTL进行富集分析

for(randomid in 1:3){
conid.random=sample(conid,150000,replace = FALSE)
con1=con[con$id %in% conid.random,]
con1=merge(con1,anno,by="id",all.x=T)

gtx=GTEx

result=data.frame(matrix(,nrow = 1,ncol =4))
names(result)=c("Term","overlap.num","OR","P.value")
result =result[-1,]


c=length(intersect(as.character(con1$id),as.character(gtx$id)))
d=150000-c
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
rt[1,]=c("control",c,NA,NA)
result=rbind(result,rt)

a=length(intersect(as.character(filea$id),as.character(gtx$id)))
b=dim(filea)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("117012",a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(filea1$id),as.character(gtx$id)))
b=dim(filea1)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("61725",a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(file$id),as.character(gtx$id)))
b=dim(file)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("53425",a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(file1$id),as.character(gtx$id)))
b=dim(file1)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("8544",a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(file2$id),as.character(gtx$id)))
b=dim(file2)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("807",a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(file3$id),as.character(gtx$id)))
b=dim(file3)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))result
rt[1,]=c("200",a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(asm1$id),as.character(gtx$id)))
b=dim(asm1)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("13649.117KASH.22OKASM",a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

rtn =paste("./GTEx.result/zong.GTEx.eQTL.",".enrich.random",randomid,".csv",sep="")
write.csv(result,rtn,quote=F,row.names=F)
}


###mQTL富集

meQTL =read.csv("/mnt/md1200/7/zhaocunyou/meQTL_by_wzj/GSE74193_jaffe_onlineTable2_meQTLs_allPairs.csv",header=T)
meQTL$id =paste(meQTL$snpChr,meQTL$snpPos,sep=":")

for(randomid in 1:3){
conid.random=sample(conid,150000,replace = FALSE)
con1=con[con$id %in% conid.random,]

result=data.frame(matrix(,nrow = 1,ncol =4))
names(result)=c("Term","overlap.num","OR","P.value")
result =result[-1,]

c=length(intersect(as.character(con1$id),as.character(meQTL$id)))
d=150000-c
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
rt[1,]=c("control",c,NA,NA)
result=rbind(result,rt)

a=length(intersect(as.character(filea$id),as.character(meQTL$id)))
b=dim(filea)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("117012",a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(filea1$id),as.character(meQTL$id)))
b=dim(filea1)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("61725",a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(file$id),as.character(meQTL$id)))
b=dim(file)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("53425",a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(file1$id),as.character(meQTL$id)))
b=dim(file1)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("8544",a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(file2$id),as.character(meQTL$id)))
b=dim(file2)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("807",a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(file3$id),as.character(meQTL$id)))
b=dim(file3)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("200",a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(asm1$id),as.character(meQTL$id)))
b=dim(asm1)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("13649.117KASH.22OKASM",a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

rtn =paste("./mQTL.result/zong.mQTL.","enrich.random",randomid,".csv",sep="")
write.csv(result,rtn,quote=F,row.names=F)
}


